#batch file for euler
#generate bootstrap distribution for Tajima D
load('grouse_mar_df3.txt')

D.CI.mat <- array(0,dim=c(4,150,10000))
n <- 50
tmp.depth.loc <- grep('depth',names(grouse.mar.df3))

for(pop in 1:4){
	for(S in 1:dim(D.CI.mat)[2]){ #no of segregating sites
		for(rep in 1:dim(D.CI.mat)[3]){
			tmp.samp <- grouse.mar.df3[,tmp.depth.loc[pop]]>=n & (grouse.mar.df3[,tmp.depth.loc[pop]-1]/grouse.mar.df3[,tmp.depth.loc[pop]])>=0.01
			tmp.data <- grouse.mar.df3[tmp.samp,c(tmp.depth.loc[pop]-1,tmp.depth.loc[pop])][sample(1:sum(tmp.samp),S),]
			tmp.data$freq <- tmp.data[,1]/tmp.data[,2]
			dummy.count <- floor(n*tmp.data$freq)
			tmp <- (n-dummy.count)*dummy.count
			k.hat <- sum(tmp)/choose(n,2)
			a1 <- sum(1/seq(1,n-1))
			a2 <- sum(1/seq(1,n-1)^2)
			b1 <- (n+1)/(3*(n-1))
			b2 <- 2*(n^2 + n + 3)/(9*n*(n-1))
			c1 <- b1 - 1/a1
			c2 <- b2 - (n+2)/(a1*n)
			e1 <- c1/a1
			e2 <- c2/(a1^2 + a2)
			D.CI.mat[pop,S,rep] <- (k.hat - S/a1)/sqrt(e1*S + e2*S*(S-1))
			if(rep%%1000 == 0) cat('population ',pop,'; S ',S,'; rep ',rep,'\n')
			}
		}
	}
rm(pop,S,rep,tmp.data,dummy.count,tmp,k.hat,a1,a2,b1,b2,c1,c2,e1,e2)
save(D.CI.mat,file='D_CI_mat.txt',ascii=TRUE)

